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Stable Evaluation of Polynomials 

C. Mesztenyi* and C. Witzgall** 

(November 3, 1966) 

A class of Newton forms 

P{x)=a Q + ai{x — x )-\- . . . -\-a n (x — x ) . . . (x — x n -i) 

are discussed which admit a stable evaluation algorithm in an interval [A, B]. Stability is defined in 
the paper. The estimate 

AP <2 + 6 M'{L)L 
\P\- z + t} M(L) ' 

where L: = B — A and M(x): = |« () | + |«i|ac+ . . . + |«n|*", is shown to hold for the relative error of 
evaluation of P{x) in [A, B]. 
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There are algorithms for polynomial evaluation 
which require fewer operations (Motzkin [4], 1 Belaga 
[1], Pan [6], Cheney [2], Knuth [3]), but these algorithms 
are prone to lose significance. The Horner scheme 
itself leaves much to be desired in this respect, and this 
fact has led to the study of methods which reduce the 
loss of significance, but may require additional opera- 
tions (Rice [7], Lawson [9]). 

An evaluation algorithm depends on the form in 
which the polynomial is given. Forms on which a 
"stable" evaluation algorithm can be based are called 
"well-conditioned." This note calls attention to a 
class of well-conditioned Newton forms. 2 These forms 
involve more parameters and require more operations 
for evaluation than the normal form described above. 
The number of operations is about the same as for 
evaluating an expansion by Tchebychef polynomials, 
but in general more parameters are required. 

Loosely speaking, an algorithm is "unstable" if 
one observes for some arguments an excessive loss 
of significance. A rigorous, but somewhat narrow, 
definition of stability will be based on an idealized 
floating point arithmetic and an error analysis similar 
to the one carried out by v. Neumann and Goldstine 
[51, however referring to the relative rather than the 
absolute error. The arithmetic will be finite in that 



it will restrict itself to numbers which are "represent- 
able" with mantissae of a fixed finite length. 3 It will 
be infinite in that it will permit arbitrarily large positive 
and negative exponents. Thus there will be infinitely 
many representable numbers. Some of the theoretical 
results of this paper hinge on this idealization, and 
might be considered artificial. Nevertheless, these 
results point in the direction of desirable improve- 
ments, and provide a reasonable classification of 
algorithms from the view point of stability. 

Our definition of condition and stability differs 
from that given by Rice [7]. 

1. Minimal Newton forms. The evaluation algo- 
rithms to be considered are based on the 



Newton form 



of polynomials: 



P(x) = ao~\- ai(x — xo)+ a2(x — xo)(x — Xi)-\- . . . 

-\-a n (x — Xq) . . . (x — Xn-\). (1.1) 



Polynomials in this form are evaluated by the follow- 
ing adaption of the Horner scheme: 
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D n : = a„ 

Di: = ai + (x — Xi)D i + 1 , 
P(x):=Do. 



,0, 



(1.2) 



3 0.53i 15 is representable in decimal arithmetic of mantissa length .s 3* 2; O.lio is not 
representable in any binary arithmetic. 



11 



Instability can arise when adding the coefficients Oi 
and when subtracting the "critical values" xu There 
are three reasons for being less concerned with the 
effect of the subtractions: 

First, if the values x and x\ are both representable 
in the given finite arithmetic, then forming the dif- 
ferences % — %i is a "stable" operation. A similar 
advantage does not obtain for the addition, even if 
the coefficients ai are representable. Furthermore, 
the effect of replacing x and x% by representable num- 
bers close to them, in other words, the effect of 
"rounding," can be estimated by means of the dif- 
ferential formula: 



dP = [Di]d(x — xo) + [D 2 (x — xo)]d(x — xi) + . . . 

. . . + [D n (x — Xo) . . .(x — Xn-2)]d{x — X n ). 



Here the D\ are the intermediate results of (1.2). 

Second, if the terms of the sum Dn = at + (x — xt)Di + i 
are of the same sign, then the effect of the relative 
error of (x — Xi)D i + 1 on the relative error of D { will 
be weighed by the ratio \{x — Xf)A + i/A|. The rela- 
tive error of x — x\ tends to be large when \x — Xf\ 
becomes small. Thus the relative error caused by 
forming the difference is toned down by the above 
ratio precisely when it tends to be high. On the other 
hand, if x and —x\ have the same sign, then one does 
not observe a similar toning down of the relative error 
incurred by the preceding addition. 

Third, we give a theoretical reason. We shall be 
able to show that x — x\ can be evaluated "stably" 
if additional coefficients and operations are intro- 
duced, provided x is representable, even if x\ is not 
representable. 

Recognizing the additions as the main cause of 
instability, the idea is to neutralize them by choosing 
the critical values x% so as to ensure that a x has always 
the same sign as {x~Xi)Di + \. 

This is indeed possible in any given interval of 
evaluation [A, B\ The first step of constructing such 
a Newton form is to divide out all zeros of the poly- 
nomial P(x) in [A, B]: 

P(x) = :(x — Xo) . • • (x — x m )P*(x). 

One then defines a =ai = . . . =a m =0, and the 
problem is clearly reduced to finding a suitable Newton 
form for P*(x). 

Let us therefore assume that P(x)t^0 in [A, B]. 
If P(x) is given in Newton form (1.1), and if a has the 
same sign as the sum of the remaining terms for all x 
in [A, B\ then sign (a ) = sign (P(x)) 9 and \a \ ^ \P(x)\ 
in [A, B\ This suggests the following definition of a () : 



Then P(x) — a has zeroes 

Xq,X\, . . . , Xk-l 

in [A, B\ Zeroes are represented as often as their 
multiplicity indicates. If, for instance, y is a double 
zero of P(x) in [A, B] then y occurs twice among the 
xu We have 

P(x) = Do(x) = ao-\-(x — x \ . . . (x — x k -i)D k (x). 

The polynomial D k (x) does not vanish in [A, B\ and 
can therefore be treated in the same way as P(x). 
The coefficient a k is denned by 

\a k \:=min{\D k (x)\ \xe[A,B]} 
sign (a k ) : = sign (D k (x)), 

and Xjc, . . . , x\-\ 

are the zeroes of D k (x) — a k . We have again 

D k (x) = a k + (x — Xk) • • • (x — xi-i)D(x), 

and so on. We call the Newton form so constructed 
the 

minimal Newton form 

of the polynomial P(x) in [A, B\ 

If all zeros of P(x) are in [A, J5], then the minimal 
Newton form of P(x) in [A , B] coincides with the "prod- 
uct form" a n Tl{x — Xi) of the polynomial. If P(x) does 
not vanish in [A, B], then all zeros of odd multiplicity 
encountered in the construction of the minimal Newton 
form are at the ends of the interval [A, B]. 

We shall show later that a minimal Newton form can 
be evaluated in a "stable" manner if the argument x is 
representable. 

2. On transformations into stable forms. It is quite 
clear that the process of determining the minimal 
Newton form is, in itself, a very unstable process, 
which has to be carried out in double or even triple 
precision. This appears to be a principle: the con- 
stants of a better conditioned form contain more 
"information," and if one wants this additional informa- 
tion, then one has to work for it. 

For this reason, we shall not concern ourselves with 
the difficulties of finding the constants of algorithms. 
These difficulties are bound to increase with the quality 
of the end-product. Besides, minimal Newton forms 
will be used only if a polynomial has to be evaluated 
repeatedly, as for instance in function subroutines. 

3. Evaluation of derivatives. The Horner scheme is 
frequently extended to determine the value of the first 
derivative of a polynomial, using the algorithm: 



|oo|: = min {\P(x)\ \xe[A,B]} 

sign(ao) : =sign(P(z)). 



E n :—D n 

Ei: = Di + xD i + u 



i = n — 1, . . . , 1, 



(3.1) 
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P'(x): = E u 

where the Di denote the intermediate results of the 
Horner scheme for- -P{x). We shall use this fact in 
section 11 for deriving an error estimate. 

Similar algorithms exist for Newton forms. If the 
D, are the intermediate results of algorithm (1.2) 
applied for a given x, then: 

P'(x) = D 1 + D 2 (x-x ) + • • • 

-\-D n (x — xo) • . . (x — Xn-2), 

which suggests the algorithm: 

E n l=D n 

Ei\ = D}-\-(x — Xi-i)E i + u i = n — 1, ... ,1 
P'(x):=E i . 

This algorithm is in general not stable, even if the 
original polynomial P is in minimal Newton form. 

4. Numerical example. 4 Consider the ill-condi- 
tioned polynomial a + . . . + .% 5 with the coefficients: 

a = 4. 10074 70239 8387 

a x =-11. 29173 84073 737 

a 2 = 8. 42475 03796 1924 

a 3 = 0. 92113 31318 58071 

a 4 = -3. 05937 81605 8204. 

Its minimal Newton form in the interval [0, 1] has the 
following coefficients and critical values: 

bo= 0.00103 19917 44066 05 

6x = 0.0 

b 2 = 3.41269 84126 9841 

A:, = -1.87912 08791 2088 

64 = 0.60784 31372 54902 

65=1.0 

*o= 0.83361 06489 18469 

Xi=X 

X2 = X\i= 1.0 
xi = 0.0. 

To compare the loss of significance incurred in 
evaluating the polynomial in normal form and minimal 
Newton form, the interval [0, 1] was divided into 50 
subintervals, in each of which 50 random arguments 
were generated. The maximum number of significant 
digits lost was plotted for each subinterval (fig." 1). 

5. Representable numbers. This and the following 
sections will contain some theoretical examinations. 
An arithmetic will be specified, with respect to which 
we shall define stability. Actual error estimates are 
then derived for the evaluation of minimal Newton 
forms. 
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Consider a floating point arithmetic with number 
base 13 and mantissa length s, and assume that there 
are no restrictions on the size of the exponents. 
Denote by 



the set of all numbers which are representable in this 
arithmetic. To each real number x, we assign a num- 
ber xeX such that no other number in £ is closer to x. 
Thus 



\x — x\ ^ \x — y\ for all jeS. 



(5.1) 



In particular, x = x if xeX. We assume that (— x) 
= — x, therefore — 2 = 2. 
With the notation 



a(x) : = 



if xeX 

1 otherwise 



one verifies that 5 



\x — x\ ^ \x\a(x)e ^ \x\e where e:=- fi s+l (5.2) 

(Scarborough's Theorem I [8] p. 3). Note that (5.2) 
holds only if the arithmetic does not restrict the size 
of the exponents. Indeed, it follows from (5.2) that: 



x = implies x = 0. 



(5.3) 



1 The computer time for this calculation was supported by the Computer Science Center 
<>f the University of Maryland under Grant NsG 398 from the National Aeronautics and 
Space Administration. 



5 Actually a sharper estimate holds: 
greatest integer smaller than log/s \x\. 



where t : = [log/i \x\\, that is. t is the 
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In other words, zero is the only number which gives 
zero when rounded. This is clearly only true if 
negative exponents of arbitrarily large modulus are 
admitted. 

6. Stable algorithms. We now assume that each 
addition, subtraction, multiplication, and division of 
two numbers in 2 is carried out exactly, and that the 
result is rounded afterwards. Furthermore, each 
given number which is not in 2 must be rounded before 
it enters an operation. These rules appear to be a 
reasonable idealization of a floating point arithmetic 
which combines two registers to hold the result of 
each operation, normalizing and rounding to single 
precision subsequently. 

An algorithm aims at computing a result rasa func- 
tion of given parameters a, 6, c, . . . , but in reality 
computes a quantity r. The error that is generated 
in the course of the algorithm then is defined by 

Ar: = \r-?\. 

Consider, for instance, the algorithm which consists 
of multiplying two given numbers u and v: 

p : = uv 



Then p: = uv, and (5.2) gives 

\p — p\= \uv — uv + uv — uv-\- uv — p\ 

^\uv\e-\-\uv\e J r\p\e. 

In view of 

M^|a| + |i*-a|^(l+€)|5|, 

|ui|^|pMttv-p|^(l + e)|p|, 

one has 

;3 + 3e + € 2 . (6.1) 



u and v: 



Ail 

«IpI 



|p— pi 



€|P| 



The upper bound of the above error term does not 
depend on u and v. Hence we say that the formation 
of the product uv from given numbers u and v is 
"stable." 

In general, we call an algorithm which computes r 
from given parameters a, b, c, . . . 

stable, 

if there exists an error estimate of the form 

Ar 



L 



er 



where k does not depend on the values of the param- 
eters a, b, c, . . . 

We turn now to the addition of two given numbers 



s = u + v. 



Then §= u~\- v, and one has after a similar calculation: 



es \s\ s\ 



(6.2) 



The estimate (6.2) can be improved. For instance, the 
first term in (6.2) can be deleted if \s\ ^ min {\u\, \v\}. 
However, no matter how refined the estimate, it will 
always depend on u and v. Indeed, u-\-v^0 and 
5 = may both hold simultaneously, precluding any 
finite bound. The addition of two numbers is therefore 
not stable. 

In the special case a(u) = a(v) = 0, that is, if both 
u and v are representable, then (6.2) proves the addition 
to be stable. This is also true if u and i;have the same 
sign. 

7. Simplifying error estimates. Henceforth we shall 
assume that e is very small. Thus we shall delete e 
in expressions of the form k + le, if both k and / are 
parameter independent: 



k+le = k. 



(7.1) 



This leads us, for instance, to equate \r\ with |r|, since 

|f|-Ae|r| = |r|(l-A*) 



r - r-r 



and similarly \r\ ^ |r|(l — ke), provided the algorithm 
producing r is stable. 

We state without proof that, if an algorithm is shown 
stable with the help of rule (7.1), then it could have 
been shown stable without it. 

8. Error propagation. Suppose that p, a, and v are 
intermediate results of some algorithm, and that 

p= uv 

according to this algorithm. We want to estimate the 
error of p in terms of Az/ and Av, and the error generated 

by rounding and multiplying. Now p = uv, and analo- 
gous to the derivation of (6.1) we have 



Ap^ |v|Att+'|ii|At;+|p|€. 



If the intermediate results u and v were arrived at in a 
stable manner, then this is also true for p. By (7.1) 
we can therefore consider \u\ and \u\, as well as \p\ and 
\p\, to be equal. Similar considerations hold for the 
division of two intermediate results. We conclude: 
If u and v are intermediate results of the algorithm 
which are obtained in a stable manner, then 



Auv ± ] 



1+ A^ + A^ 
e\u\ e\v\ 



(8.1) 
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The estimate (8.1) is usually postulated without 
assumptions about the algorithm. It is argued that 

An Av 

l — r and 7-7 are comparable to e in magnitude, ... a 

\u\ t \v\ 

claim that comes close to assuming stability of 
the algorithm . . . , and that therefore the term 

An Av . , . 

— ; — r • — j— 7 € may be neglected. 

e\u\ e\v\ 

The following estimate can be derived without 
stability assumptions: 



A(u±v) 



1+- 



e \u : 



AtL + A 
e\u\ \u: 



Av 
e\v\' 



(8.2) 



9. Stability of minimal Newton forms. Let u» now 
turn to algorithms for the evaluation of polynomials. 
Once and for all we will assume that 

(9.1) the argument x for which a polynomial P(x) is to 
he evaluated is representable. 



since x is representable. Hence 

^f** (9.2) 

€\x-y\ 

The algorithms ({x — y)-\-s) — r and (."?— 7) + (s — r), 
where x is not required to be representable and where s 
denotes the remainder term x — x of X, are not stable. 
Indeed, if x = y, then both algorithms consist of simply 
subtracting r from s, which was seen to be unstable in 
section 6. 

Consider the algorithm (1.2) for evaluating Newton 
forms. Since all other operations of (1.2) are stable, 
the entire algorithm (1.2) is stable provided the dif- 
ferences x — Xi are evaluated in the stable manner 
described above. 

For the step 

Di : = ai J r(x — Xi)D i + i, 
we have by (9.2), (8.1), and (8.2) 



Consider the evaluation of a polynomial in minimal 
Newton form by the adapted Horner scheme (1.2). 
All operations in this algorithm are stable except the 
subtractions x — xi. The existence of a stable evalua- 
tion thus depends on the possibility of evaluating the 
differences x — X\ in a stable manner. 

To evaluate the difference x — y, we split y into its 
representable part and remainder 



y = y+r, 

and evaluate (x — y)—r. Since i_and y are both 
representable, and therefore A.r = Ay=0, we have by 

(8.2): 



A(jr-y) „ 

e\x-y\ " 

Then again by (8.2): 

A(x-y) ^ \x-y\ A(x-y) 



1. 



e\x-y\ 



\x-y\ e\x-y\ 



A/ 



1 + 



l*-y| + H 



x-y c /■ 



In view of |# — y| = \x — y-\- r\ ^ \x — y\ + \r\ this gives 



A(x — Xj)Dj + i 
e\(x — Xi)D i + 1 



1 A(x-xj) AD i+1 c5 AD i+] 



ex 



e\D i + ] 



e\Du 



M)i \at\ Am \(x-Xi)I) i+ , 1 A(x-Xi)D i + i 

elAP \Di\'e\a,\ \Di\ " e\(x-Xi)D i+l \ 



< \a,\ \(x- Xi )D i + i \ I 1 Agj+J 

" IAI IAI \ e\D i+] 



\ai\ + \(x-xt)Dt + i\ , \(x-xt)D i + i\ ( A , M) i + 



V e|A + i|/ 



Since minimal Newton forms are constructed so as to 



Ne 
ensure that |«,-| + \{x — xt)Dt + i\ = |Di|, we finally have 



AD,- _ 2 + l(.y-.v,-)A,- + il ^ + AD, 



«|A| ~~ ' |A| V ' 6|D»+i| 

from which one immediately concludes that 



AD, _. , AD» + , 

+ 



(9.3) 



€\D,\ 



e\D t+ i\ 



Since AD„ «£ |D„|e, we have finally 



A(*~y) <2| 2|r| 
e|*-y| ~~ i-v — >-| 



AP 

€|P| 



6n + l, 



(9.4) 



Now according to (5.1) 



■y\, 



where n is the degree of the polynomial P(x). 

Note that application of the estimate (8.2) is justi- 
fied since the factors of the products (x — Xi)Di + i have 
been arrived at in a stable manner. 
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10. Product forms. Every polynomial with real co- 
efficients can be written in product form; 



P(x) = a n l\(x-Xi) J] {di + {x-XiY), 



k+l-l 

n 

i=k 



di>0, fc + 21-- 



(10.1) 



A stable evaluation algorithm results if the differences 
x— xi are evaluated in the stable manner described 
in section 9. The algorithm involves not more than 
2/i+l parameters, 2n additions, and n multiplications. 
The following error estimate follows from (9.2), (8.1), 
and (8.2): 



AP 

e|P| 



5A -h 11/+ 1 ^6n. 



(10.2) 



11. Improved estimate. The reader observes not 
only, that the error estimate (10.2) is better than the 
estimate (9.4) for minimal Newton forms, but also that 
fewer operations are required. It is therefore ques- 
tionable whether minimal Newton forms should be 
considered at all, unless the error estimate (9.4) can 
be improved. We proceed to show that this is indeed 
possible. 

Let L denote the length of the interval [A, B] and 
define 

A:=H+|oi+i|£+. . .+|a n |I>-\ i = 0, . . ., n. 

Then \D t \^Bi in [A, B] since 

Di = ai + (n+i(x — Xi)-\-. . .-\~a n (x — Xi) . . . (x — x n -i). 



Therefore 

\(x-Xi)D i+i \ 



Di 



and (9.3) reduces to 



I A- 1 



I".- 1 
B, 



B, 



ADi -2 + ^U + ^)L. 



eA 



B, \ e\D l + l \ 



Consequently 



e\Pr 2+ B L+ BSi 



and therefore 



,+ 



B\ . . . B n 



Bo 



• B n -\ 



D\ 



€\P\ 



2 + 7 -(^ 1 L + J e 2 L 2 +. 



.+B n L"). (11.1) 



If we introduce the polynomial 

M(L):=|a |+|ai|L+|a 2 |L 2 +. . . + \a n \L n , 

Then the Bi are the intermediate results of the Horner 
scheme applied to M(L). Thus B = M(L). More- 
over, we obtain from (3.1) that B x + . . . +B„L n ~ 1 
= M'(L), and we may rewrite (11.1) as follows: 



AP 

e|P| 



2 + 6 



M'(L)L 

M(L) ' 



(11.2) 



This error estimate establishes a preference for 
minimal Newton forms in small intervals. 

Polynomials with nonnegative coefficients are 
minimal Newton forms for all intervals [0, x] with 
x > 0. For such polynomials one derives analogously 



AP 

e|P| 



2 + 2 



P'(L)L 
P(L) ' 



12. Approximately minimal Newton forms. Each 
sequence of n critical values %\ determines a unique 
Newton representation of a given polynomial P(x). 
If the values xf approximate the points %\ which deter- 
mine the minimal Newton form, then the Newton form 

P(*) = a* + a*(x -**)+. . . 

+ «*(*-**) . . . (*-**_!), 

which is determined by the values xf is approximately 
minimal, and for all practical purposes, it is stable. 
For certain polynomials of low degree approximately 
minimal Newton forms may be even rigorously stable 
in the sense of section 6. Consider for instance the 
quadratic polynomial in minimal Newton form 



Q(x) = «o+ (x — z) 



a > 0. 



The Newton representation of Q(x) determined by the 
rounded value z takes the form: 



Q(x) = b + b 1 (x-z)+ (x-z) 2 . 



(12.1) 



where bo>r 2 , b\= — 2r, with r:=z — z. We proceed 
to prove that the evaluation of (12.1) by algorithm 
(1.2) is stable. We put s :=x — z. Then 

D:=b l +(x-z)=-2r + s 
Q =b +{x-z)D=b Q + sD. 
We have 

AD^2|r|e + s|e|+|Z)|e 
AsD ^ 2\sD\e+\s\AD ^ 3\sD\e + 2\rs\e + s 2 e 



16 



and 



Hence 



AQ^\b \e + 3\sD\€ + 2\rs\e + s 2 6 + Qe (12.2) 
By (5.1), \s-r\=\x-z\^\z-z\=\r\. Hence 



Q = ao+ (x— -a 



(b -r*) + (s-r) 2 ^b 



for all representable x. This in turn implies 
Q 22 sD ^ 0. As a consequence, we may delete most 
of the absolute bars in (12.2): 

— ^bo + SsD + 2|rs|+s 2 + Q^2Q + 2sD + 2\rs\+s 2 . 



If rs ^ 0, then 2|rs| + s 2 = s 2 — 2rs = sD. Hence 



— ^2Q + 3sD^5Q. 

e 



If rs > 0, then |s|^2|r|. Indeed (s — r) 2 ^ r 2 implies 
2rs ^ s 2 , and if rs > 0, then we have also 2|r||s|^|s| 2 . 
Furthermore, sD + 2\rs\= sD + 2rs = s 2 . Finally, we 
note that 



Q (s-r) 



- 2 ^4for\s\^\2r\. 



e|<?f 



2 + ^> + 



2* 2 



Q (s-r) 2 



11. 



This estimate is as good as the one derived from (10.2) 
taking into account the fact that a 2 = 1. 
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